Part of the Accelerator Readiness Review is looking at muons that leave the ring. To this end, we (Renee, Nathan, James, and myself) have a preliminary version of MDC-2 with full injection, unified fields, and “ghost” detectors in a Geant parallel world (so these sensitive volumes can be placed without overlapping physical structures). We can remove various magnets and fields and see how that affects muons escaping the beam.
# Load libraries
library(readGallery)
library(stringr)
library(dplyr)
library(ggplot2)
library(rgl)
I have run preliminary MDC-2 (using Renee’s FCL files), making 10K events for various scenarios. Files are currently stored in /pnfs/GM2/scratch/users/lyon/arr_20170313 and directories within. Ghost detectors included are Renee’s cylinder just at the outer vacuum wall and encompossing the inflector as well as my ghost detector that is a cylindrical shell just on the inside of the world cube. Note that this code includes the magnet yoke steel.
# Let's look at only the 180* job series (so we just get 10,000 events per scenario)
system('ifdh ls /pnfs/GM2/scratch/users/lyon/arr_20170313/180*/*.root | grep .root', intern=T) %>%
xrootify() -> arrFiles
arrFiles
[1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"
[3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"
[4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"
[5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"
[6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"
Let’s pluck out the scenario.
scenarios <- str_match(arrFiles, 'arr_20170313/.+_unified_(.+)_cyl')[,2]
scenarios
[1] "noKickerQuads" "noKicker" "noQuads" "noInflector" "everything" "noDipole"
What data are in these file? Well, product_sizes_dumper doesn’t seem to work across XRootD. That’s too bad.
Here’s a picture from Root.
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
Create a reader class (modify outside of this notebook) for TrackingActionArtRecord
readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')
trackingActionReaderFor documentation, here is the TrackingActionArtRecordReader reader class
readr::read_file('trackingActionReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag, evLimit):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz']
self.evLimit = evLimit
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
if ev.eventEntry() > self.evLimit:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Let's only accept the muon
if e.trackID != 1 or e.parentTrackID != 0:
continue
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID, e.volumeUID,
e.status, e.x, e.y, e.z, e.px, e.py, e.pz ])
return True
Make an instance
trackingActionReaderClass <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader
Load the data - let’s just look at one file since it takes a long time
taReader <- trackingActionReaderClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[1], taReader)
Opening first file...
'cannot obtain credentials for protocol: Secgsi: ErrParseBuffer: error getting user proxies: kXGS_init: unable to get protocol object.'
Error in py_call_impl(callable, dots$args, dots$keywords) :
Exception: gallery::Event::gallery::Event(const vector<string>& fileNames, bool useTTreeCache = true, unsigned int eventsToLearnUsedBranches = 7) =>
---- FileOpenError BEGIN
Failed opening file 'root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root'
---- FileOpenError END
(C++ exception of type cet::coded_exception<art::errors::ErrorCodes,&art::ExceptionDetail::translate>)
Detailed traceback:
File "/home/me/R/x86_64-redhat-linux-gnu-library/3.3/readGallery/python/readGallery.py", line 88, in getGalleryData
ev = ROOT.gallery.Event(fileVector)
As of right now, the dCache XRootD door is done.
nrow(taDF)
[1] 120000
Set the scenario (we’ll need that again, so make it a function)
setScenario <- function(df) {
df %>% mutate(scenario = factor(fileEntry, levels=0:(length(scenarios)-1), labels=scenarios))
}
taDF %>% setScenario() %>% select(scenario, everything()) -> taDF
taDF
How many track action hits do we get per scenario.
taDF %>% group_by(scenario) %>% tally()
Let’s just get the death of the muons.
taDF %>% filter(status == 1) -> taDeathDF
We need the volume ID. I have a art FCL to run to get this information. I have files and they each may have a different set of volume IDs. Let’s try to do this in parallel. The examples in the help for parallel::mcparallel seem to be useful here.
Create the command strings,
runForVolIDString <- function(i) {
str_interp(
"PVS_CSVOUT=${csvout}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
list( csvout=scenarios[i], fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrFiles[i])
)
}
runForVolIDStrings <- sapply(1:length(arrFiles), runForVolIDString)
runForVolIDStrings
[1] "PVS_CSVOUT=noKickerQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[2] "PVS_CSVOUT=noKicker_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"
[3] "PVS_CSVOUT=noQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"
[4] "PVS_CSVOUT=noInflector_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"
[5] "PVS_CSVOUT=everything_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"
[6] "PVS_CSVOUT=noDipole_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"
Let’s run in parallel
jobs <- lapply(1:length(arrFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)
Now we need to load the csv files. Let’s do that in parallel too
jobs <- lapply(1:length(arrFiles), function(i) paste0(scenarios[i], "_volNames.csv") %>% read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
volNameTables <- mccollect(jobs)
names(volNameTables) <- scenarios
Let’s see what we got
volNameTables[["everything"]]
Join them together when the scenario name.
volNameTable <- map2_df(volNameTables, scenarios, function(df, sce) df %>% mutate(scenario = sce) )
volNameTable
Now let’s do lookups… [Note that I can’t do what I originally thought, which was split up the data frame into groups by scenario and then put in the factor by group – because in the end the rows get rejoined and then row-binded factor gets messed up.]
taDeathDF -> hold # in case we mess up
taDeathDF %>% inner_join(volNameTable) %>% select(scenario, eventEntry, volName, everything()) -> taDeathDF
Joining, by = c("scenario", "volumeUID")
joining factor and character vector, coercing into character vector
taDeathDF
We can check this – we should see different volume names associated with volume IDs for different scenarios
taDeathDF %>% distinct(scenario, volumeUID, volName) %>% arrange(volumeUID)
Let’s plot where things die
taDeathDF %>% group_by(scenario, volName) %>% tally() %>% arrange(scenario, desc(n)) -> deathVolumeTable
deathVolumeTable
Did any of them actually decay?
taDeathDF %>% filter(volName == "ArcSection[00]")
Uh oh!!!! Hardly any decay!
taDeathDF %>% ggplot( aes(x = volName) ) + geom_bar() + theme(axis.text.x = element_text(angle=90, hjust=1))
Let’s look at the ‘everything’ scenario.
everythingFile <- which(scenarios == "everything")
I made a new TrackingAction reader
readr::read_file('trackingActionReaderCheckDecay.py') %>% cat
# Store the muon at death and whatever it decayed to at birth
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionWithDecayReader(GalleryReaderBase):
def __init__(self, inputTag, evLimit):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz', 'e', 'turn']
self.evLimit = evLimit
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
if ev.eventEntry() > self.evLimit:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Let's only accept the muon at death or whatever the muon decays to at birth
if (e.trackID == 1) or (e.parentTrackID == 1 and e.status == 0):
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID, e.volumeUID,
e.status, e.x, e.y, e.z, e.px, e.py, e.pz, e.e, e.turn ])
return True
Make an instance
taDecayClass <- createReaderClass_from_file('trackingActionReaderCheckDecay.py')$TrackingActionWithDecayReader
Load the data - let’s just look at one file since it takes a long time
taReader <- taDecayClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[5], taReader)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root
tadDF <- galleryReader_df(taReader)
nrow(tadDF)
[1] 52842
Make things easier to look at
tadDF %>% mutate(r = sqrt(x*x+y*y)) %>% select(trackType, volumeUID, status, r, e, turn) -> tadDF
So I see an electron in alomst every case. But the muon appears to stop in the iron. Let’s just plot x,z
taDeathDF %>% ggplot(aes(x=z, y=x)) +
geom_point() + xlim(c(-10000, 10000))
tmp <- open3d()
ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
radius=7112,
sides=20, closed = F)
clear3d()
with(taDeathDF,
plot3d(x=x, y=y, z=z, type='p', main='Muon death positions',
xlab="x", ylab="y", zlab="z", xlim=c(-9000, 9000), zlim=c(-9000, 9000)))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
This really doesn’t look like muons are being stored.
useDataProduct('std::vector<gm2truth::GhostDetectorArtRecord>')
ghostCylinderInputTag = artInputTag('artg4:GhostCylinderDetector')
Create a reader class (modify outside of this notebook) for GhostDetectorArtRecord
readerClassSkel('gm2truth::GhostDetectorArtRecord', writeFile = 'ghostDetectorReader.py')
Here’s the reader class
readr::read_file('ghostDetectorReader.py') %>% cat
Make an instance
ghostDetectorReaderClass <- createReaderClass_from_file('ghostDetectorReader.py')$GhostDetectorArtRecordReader
ghostDetectorReader = ghostDetectorReaderClass(ghostCylinderInputTag)
Load the data
getGalleryData(arrFiles, ghostDetectorReader)
gcDF <- galleryReader_df(ghostDetectorReader)
How many events did we get?
nrow(gcDF)
head(gcDF)
Make factors…
Some common particles are
# See http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf for PDG value meanings
pdgs <- c('e-'= 11, 'nu_e'= 12, 'mu-'= 13, 'nu_mu'= 14, 'tau-'= 15, 'nu_tau'= 16,
'e+'=-11, 'anti_nu_e'=-12, 'mu+'=-13, 'anti_nu_mu'=-14, 'tau+'=-15, 'anti_nu_tau'=-16,
'gam'=22, 'p'=2212, 'n'=2112, 'anti-p'=-2212, 'anti-n'=-2112)
pdgs
gcDF %>% mutate(scenario = factor(fileEntry, levels=0:(length(scenarios)-1), labels=scenarios),
pdg = factor(particleID, levels=pdgs, labels=names(pdgs))) %>%
select(scenario, eventEntry, pdg, everything())-> gcDF
gcDF
How many lost per scenario
gcDF %>% filter(pdg == 'mu+') %>% group_by(scenario) %>% tally()
How many events?
gcDF %>% group_by(scenario) %>% tally()
We can see the muons that escape the vacuum chamber by looking at the ghost detector Hits from Renee’s detector (has art input tag of artg4:GhostCylinderDetector). We can also look at hits from my ghost detector that sits just inside the world cube.
Tell PyRoot to load the necessary dictionaries.
readGalleryPy$provide_get_valid_handle('std::vector<gm2truth::GhostDetectorArtRecord>')
gh_cyl_tag <- readGalleryPy$createInputTag('artg4:GhostCylinderDetector')
gh_nwd_tag <- readGalleryPy$createInputTag('artg4:GhostNearWorldDetector')
Let’s try to load data!
pyClass = '
class GhostDetectorArtRecordReader:
"""Reader object to read Ghost Detector Art Record objects"""
def __init__(self, inputTag):
self.vals = []
self.inputTag = inputTag
self.getValidHandle = None
self.names = ["fileEntry", "eventEntry", "particleID", "trackID", "parentTrackID", "x", "y", "z",
"px", "py", "pz"]
def colnames(self):
return self.names
def prepare(self, ROOT, ev):
self.vals = [] # Protect againt re-run
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.GhostDetectorArtRecord))
def fill(self, ROOT, ev):
gh_cyl_h = self.getValidHandle(self.inputTag) # Get the valid handle
if not gh_cyl_h.empty(): # Does it have data?
gh_cyl = gh_cyl_h.product() # Get the corresponding data product vector
for g in gh_cyl: # Loop over elements and fill
if g.trackID == 1 and g.parentTrackID == 0:
self.vals.append(
[ev.fileEntry(), ev.eventEntry(), g.particleID, g.trackID, g.parentTrackID,
g.position.x(), g.position.y(), g.position.z(),
g.momentum.x(), g.momentum.y(), g.momentum.z()
])
return True
'
py_run_string(pyClass)
ghCReader <- main$GhostDetectorArtRecordReader(gh_cyl_tag)
ghNReader <- main$GhostDetectorArtRecordReader(gh_nwd_tag)
readers <- readGalleryPy$GalleryReaders(list(ghCReader, ghNReader))
Read the data
readGalleryPy$getGalleryData(fileVec1, readers) -> times
Check out the time
times$allTime
ghcdf <- as.data.frame(do.call(rbind, ghCReader$vals))
ghcdf[] <- lapply(ghcdf, unlist)
names(ghcdf) <- ghCReader$colnames()
ghcdf
format(nrow(ghcdf), big.mark = ',')
ghndf <- as.data.frame(do.call(rbind, ghNReader$vals))
ghndf[] <- lapply(ghndf, unlist)
names(ghndf) <- ghNReader$colnames()
ghndf
ghCReader$colnames()